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Abstract 


A dynamic PEM fuel cell model has been developed, taking into account spatial dependencies of voltage, current, material flows, and temperatures. 
The voltage, current, and therefore, the efficiency are dependent on the temperature and other variables, which can be optimized on the fly to achieve 
optimal efficiency. In this paper, we demonstrate that a model predictive controller, relying on a reduced-order approximation of the dynamic PEM 
fuel cell model can satisfy setpoint changes in the power demand, while at the same time, minimize fuel consumption to maximize the efficiency. 
The main conclusion of the paper is that by appropriate formulation of the objective function, reliable optimization of the performance of a PEM 
fuel cell can be performed in which the main tunable parameter is the prediction and control horizons, V and U, respectively. We have demonstrated 
that increased fuel efficiency can be obtained at the expense of slower responses, by increasing the values of these parameters. 


© 2007 Elsevier B.V. All rights reserved. 
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1. Introduction 


Fuel cell models of various levels of complexity have been 
suggested, describing their performance under an array of con- 
ditions [1-3]. These models have been used to evaluate optimal 
schemes of external heating, water management and fuel com- 
position. The regulation of the transient response of fuel cells is 
important for vehicular applications, since the power demands 
fluctuate, and the fuel cell will not be working at the optimal 
steady state designed by its manufacturer. It is desirable to con- 
trol the fuel cell so that acceptable response time for the power 
demand is ensured, while achieving high efficiencies. Dynamic 
models facilitate the design and testing of control systems. To 
this end, a dynamic empirical model for the transient response 
of a fuel cell was developed by Amphlett et al. [4]. This is 
a lumped model with no spatial dependence. Kang et al. [5] 
present an analysis of a dynamic model for a molten-carbonate 
fuel cell (MCFC) where the system is modeled as a collection 
of first order transfer functions with dead times. More recently, 
load distribution control in hybrid vehicles has been investigated 
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[6,7], motivated by the benefits in improved dynamic response 
of systems combining fuel cells and a battery. In addition, con- 
trol of the power conditioning unit is addressed [8] along with 
control of the fuel delivery system [9]. These studies typically 
treat the fuel cell itself as a type of disturbance and its control 
is not addressed. Lauzze and Chmielewski [10] attempt to per- 
form multivariable control using separate linear PID controllers. 
However, the system’s high degree of interaction between the 
different control loops is evident. 

In previous work [11], we described a framework to control 
the fuel cell using nonlinear model-based control. We pro- 
posed a transient model, based on the ideas presented by Yi 
and Nguyen [2], modeling a fuel cell along its channel. Both 
models account for heat transfer between the solid and the two 
gas channels, and between the solid and cooling water. In addi- 
tion, the water content is modeled, accounting for condensation 
and evaporation, water drag through the membrane, and water 
generation at the cathode. However, we account for the tran- 
sients in the energy balance on the solid, with all the other 
equations assumed to be at quasi-steady-state in equilibrium 
with a given solid temperature profile. This spatially dependent 
model is referred to as the “full-order” model, and is used to rep- 
resent the true process in closed-loop simulations. Furthermore, 
a first-order, time dependent model of a fuel cell has been devel- 
oped, which is fast enough to use for on-the-fly optimization of 
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solid-gas heat transfer area per unit length along 
channel (cm) 

matrix for linear constraints in NLP 
solid-coolant heat transfer area per unit length 
along channel (cm) 

vector for linear constraints in NLP 

target function weight matrix of NLP 

solid heat capacity (J (g Ky) 

diffusion coefficient (cm? s7!) 

area of current per unit length along channel (cm) 
cross-section of solid in direction of reactant flow 
(cm”) 

Faraday constant (C s7!) 

channel width (cm) 

enthalpy of overall reaction (J mol~') 

enthalpy of water condensation (J mol~!) 
current density (A cm~?) 

exchange current density (A cm~?) 

objective function for MPC 

channel length (cm) 

vector of optimization variables 

molar flow (mol s~!) 

power density (W cm~?) 

partial pressure of species j (atm) 

gas constant (8.314 J (mol K)“!) 

weight coefficient matrix for control moves 
membrane thickness (cm) 

temperature (K) 

control variable 

change in control variable value 

heat transfer coefficient (W (cm? K)~!) 
prediction horizon 

cell voltage (V) 

open circuit voltage (V) 

weight coefficient matrix for setpoint tracking 
output variable 


Greek letters 


ô length of diffusion layer (cm) 

y minimum value for gradient 

n efficiency 

Qg minimum stoichiometric ratio between hydrogen 
and current 

A waste value (mol s~!) 

p solid density (g cm7?) 

o membrane conductivity (Q cm)! 

Weft weight coefficient matrix of efficiency 

Subscripts 

a anode 

avg average 

c cathode 

cool coolant 

8 gas 


H2 hydrogen 


i time step 

in inlet 

inf surroundings 
N2 nitrogen 

O2 oxygen 

oc open circuit 
s solid 

set set point 

w water 
Superscripts 

l liquid 

sat saturation 

v vapor 


operating parameters to ensure convergence to required power. 
This reduced model was used by a model-predictive-controller 
(described below) to regulate the power output of the fuel cell. 
As demonstrated in Golbert and Lewin [11], model predictive 
control (MPC) relying on this model is more robust than stan- 
dard linear control, especially in regions of high power density. 
In addition, multiple degrees of freedom can be used to improve 
performance. 

The objective of this paper is to demonstrate that MPC can 
be further exploited to achieve both robust performance as well 
as improved fuel efficiency. We begin by reviewing the model 
predictive control method and explain the solution procedure. 
Next, we review the reduced-order fuel cell model utilized 
by the controller. Using this reduced model we describe the 
MPC formulation for fuel efficiency. Finally, we present the 
results obtained, comparing the performance with and with- 
out accounting for optimal efficiency. The results show that 
advanced model-based control can improve efficiency by uti- 
lizing the degrees of freedom in the fuel cell operation. A novel 
definition of the MPC target function is presented which simul- 
taneously achieves convergence and maximizes efficiency as 
opposed to a trade-off between the two concerns. 


2. Modeling and analysis 
2.1. Description of MPC 


Model predictive control [12] is part of a family of 
optimization-based control methods, which are based on on-line 
optimization of future control moves. Using a process model, 
the optimizer predicts the effect of past inputs on future outputs. 
Then, using the same model, it computes a sequence of future 
control moves to minimize an objective function that including 
penalties on the trajectory of predicted tracking error. The first of 
the future control moves is implemented, and the entire optimiza- 
tion is repeated from the next step on, and so on, ad infinitum. 
Feedback is used to account for the model’s inaccuracies and to 
ensure convergence. 
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The work described in this paper solves the optimization 
problem implicit in MPC simultaneously using the approach of 
Biegler [13], briefly reviewed next. A typical nonlinear discrete 
system is defined as 


Veg = LO Lo) (1) 


To simplify nomenclature, we assume no system memory, 
that is, yz; depends explicitly only on the system values at time 
k (i.e., no time lag), noting that dependencies on historical data 
can be included with no loss of generality. For such a system, a 
nonlinear problem is defined as 


min J(Au;,..-, Auy Ypo Yy) 
Auy,.--, Auy 
Yr Dy 
v U 
=X y Wo + Y Aus Au; (2) 
i m 
i=l i=l 
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Unin S Uy = Umax? Almin s Au, = AU max 


where the value of the control variables is simply the sum initial 
values and the changes in the variables: 


k 
Uk = uy + > Au; (3) 
i 


The target function, J, expresses the trade-off between the 
convergence to a desired trajectory in the outputs and the 
required moves of the control variables, u. As will be shown 
later in this paper, the target function determines the closed loop 
performance of the system. 

In reality, the control variables (throughout the control 
horizon) are the only independent variables, whereas the opti- 
mization variables are defined as the control variables, Au, and 
the output variables, y. The apparent discrepancy in the num- 
ber of degrees of freedom is resolved by the nonlinear equality 
constraints. These constraints ensure that the output values of 
the optimization solution are feasible as defined by the system 
equations given as Eq. (1). The equality constraints reduce the 
degrees of freedom of the optimization to be identical to the 
number of control moves. 

To achieve the classic form of an NLP, the optimization 
variable, m, is defined as the concatenation of the changes 
in the control variables (over the control horizon) and the 
values of the output variables (over the prediction horizon): 
m = [Auj,..., AUU, y1,.---, yv]'. The replacement of the lin- 
ear bounds in Eq. (2) by linear constraints and bounds on m is 
straightforward resulting in the final NLP which is solved at 
every time step: 


minJ(m) = m'Cm 
Pe (4) 


st. gim)=0, Am<B, Mpmin <M < Mma 


2.2. Reduced-order model 


The model used in the MPC framework is a time dependant, 
lumped parameter model of the flow channels, and membrane, 


including temperatures and water content. Since this model has 
been detailed in Golbert and Lewin [11], we state only the key 
equations here. Assuming quasi-steady-state for most of the state 
variables leaves the temperature of the solid as the only dynamic 
variable: 


dT; 
dt = fr.(Ts, Tools -+> Iavg) 
1 Usa U, 1b 
= (HE +r- 2T;) F ee Tool = Ts) 
PsC ps f f 
ee na teed An (Ts) 
f OF cell avg FL vapl £s 
l 1 1 l 
x (Mya My ain a My a M\y cin) 
QU ing 
= = (Ts = TaD) (5) 


The Nernst equation defines the dependence of the voltage 
on the current density accounting for overpotentials: 
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where the concentration overpotential is accounted for by defin- 
ing: 
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The solution of Eq. (6) gives the current density. The con- 
centrations of the reactant hydrogen and oxygen and the water 
product are governed by simple mass balances where their pro- 
duction/consumption rates are dependent on the current density. 
The spatially dependent model is simplified to enable rapid cal- 
culation for control and optimization purposes by lumping the 
spatial dependencies, which results in simple algebraic equa- 
tions (see Golbert and Lewin [11], for details): 


hL 
Mp, = OF avg F My, in (7) 
AL 
Mo, _ AF avg F Mo,,in (8) 
U,al(T, — Ty) 
Tk = Thin + -En 4, k=a,c 9 
k = Tk,in SMC, M; (9) 
Ucoolb L(Ts — Teoot) 
Teool = ‘£cool,in +F = à ve (10) 
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Since the system is assumed to be at quasi-steady-state, the 
water vapor in each channel is assumed to be at equilibrium 
with liquid water if present (i.e., at the saturation pressure), in 
which case the amount of liquid is a balance of water entering 
and exiting the channel: 


C2 
wa Pal Psat(Ta) — 1 


l Şi hL 
+ M. -Mya -eT 


w,a,in 


1 
My. = MY 


w,a,in 


(11a) 


If this results in a negative value for the liquid water, then the 
anode is not in equilibrium, and there is no liquid water present: 


+M! 


w,a,in 


v v hL 1 

Myw, a = Me yin ad Mwa = 9 (11b) 
The water content at the cathode is similarly calculated 

including the expression for water generation from the reaction. 


If the cathode water is at equilibrium: 


v Mo, + Mn, 
Hy Po/ Psat(Ta) — 1’ 


Mo = Ment Moin Mere Ir (12a) 
al a : F 2F 
Otherwise there is no liquid water present: 
M}, o atenta es "i 
l na i F 2F 
My. =0 (12b) 


In Egs. (11a), (1 1b), (12a) and (12b), a, the water drag ratio, is 
a function of the temperature, pressures and water content [14]. 
Since the water content is, in turn, a function of a Eqs. (11a), 
(11b), (12a) and (12b) define an implicit equation, the solution 
of which results in the value of a: 


a= fol(Ta, Te, My, (a), M}, (œ), Pa, Po) 
= falTa, Te, œ, Pa, Po) (13) 


2.3. Definition of optimization problem 


In this paper we will address a system where the control vari- 
ables are the dry hydrogen flow rate (assumed to be humidified), 
the coolant temperature and the average current density. Origi- 
nally, the system was defined using the solid temperature, cell 
voltage and water drag ratio, a: 


Mp,,i Tsi 
u; = | Teool,i | > y= Pi (14) 
Tayg,i Qi 


2.3.1. Objective function 

In Golbert and Lewin [11] we demonstrated that nonlinear 
MPC can satisfy changes in load demands robustly. We showed 
that the usage of multiple manipulated variables can improve 
the response time of the system, with the target function to be 


minimized being the sum of square errors from the setpoint, 
with a penalty on the moves required in the manipulated vari- 
ables. Clearly, there is potential for the optimizer to exploit the 
degrees of freedom inherent in the fuel cell design to improve the 
fuel-efficiency. In this regard, efficiency is defined as the ratio 
between the actual power produced and the heat of formation of 
the water produced if all the hydrogen feed is consumed: 


hLP 


7 15 
AHMu, ts) 


n 

Since the MPC solves a minimization problem, a waste vari- 
able is defined, 1 — n, which is to be minimized by the optimizer. 
The most obvious way to improve the efficiency is to lower the 
feed flow rate. This is offset with the need for sufficient hydrogen 
concentration to achieve satisfactory voltage. This is the rea- 
son that the controller will not reduce the hydrogen too much, 
since the concentration overpotential will become unbearable 
and compromise the power output. The optimization problem is 
now defined as the weighted sum of the performance (the dif- 
ference between the set point and the actual power), the size of 
the control steps and the local efficiency over time: 


J(Au, y) = (1 — oP — Pset) W (P — Pset) 


+o — n)"(1 — n) + Au’ SAu (16) 


where the power density, P, is not calculated, but merely is a 
component of vector y. The parameter weft expresses the desired 
trade-off between performance and efficiency, while the coeffi- 
cients of the diagonal matrix W_ are selected to scale the terms 
in proportion to their importance. 

Unfortunately, the definition of efficiency in Eq. (15) can lead 
to numerical problems, since it involves division by the hydrogen 
flow rate, which at low values, will lead to excessively large 
objective function gradients, and at high values to gradients in 
the efficiency contribution to J that approach zero. Consequently, 
either the optimizer will give excessive attention to efficiency, 
or will completely ignore it, depending on the local value of 
the hydrogen flow rate. To avoid these problems, the following 
expression for the efficiency loss, A, is used for the optimization: 
A=M n P 17 
11 — MAD AH- ( ) 

The second term in Eq. (17) is the minimum hydrogen nec- 
essary for the actual amount of power produced based on the 
change in enthalpy due to the oxidation of hydrogen. Note that 
if the hydrogen flow rate exactly matches the amount needed 
to satisfy the power demand, A = 0, and the influences of the 
power and hydrogen feed are retained in this formulation of the 
system inefficiency. Using this modified loss efficiency term, the 
objective function is now: 


J(Au, y) = (1 — oP — Ps) W (P — Pset) 


+ wet AW, A + Au’ SAu (18) 


In Eq. (18) it is noted that the coefficients of the diagonal 
matrix W, are set to 10! to scale the second term with the 
others. 
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2.3.2. Equality constraints 

The equality constraints, representing the discrete system, 
include the difference equations for the solid temperature based 
on Eq. (5) for each of the V time steps in the prediction horizon: 


8i, Y) = Ts,i+1 — Ts,i — fr,(Ts,i, Veet,i, Hi, .- JAt= 0, 
i=0,1...,V—1 (19) 


Note that the numerical solution of Eq. (5) is implemented 
using the explicit Euler method, although increased accuracy is 
possible by invoking higher-order methods. Note that for sim- 
plicity, all variables assumed to be constant (for example the 
channel pressures and inlet temperatures) are not noted in the 
equations. The equality constraint for the power density requires 
that the assumed value for the power equals the calculated value, 
where the voltage, Vcen,;, is calculated using Eq. (6): 


82(ui, Yi) = Pi — Li Veen,i = 9, i=0,1...,V—-1 (20) 


The nonlinear equation (13) for œ for each time step is also 
represented by a constraint: 


83(Ui, Yi) = Qi — falTs,i, æi, ...) = 0 (21) 


Itis well known that the fuel cell power—current curve exhibits 
a maximum value, while the voltage—current curve is nonlinear 
but monotonically decreasing. As a result of these characteris- 
tics, it is desirable to limit the operation region of fuel cells to 
always lie to the left of the peak in power, since this will ensure 
operating at high voltage and higher efficiencies. Unfortunately, 
at high power demands, the system can cross the maximum, and 
as was shown by Golbert and Lewin [11], model inaccuracies can 
lead to mistaken predictions of the sign of the gain between the 
power and current density. A controller relying on this prediction 
would lose the ability to control the system. On the other hand, 
using the voltage as a state variable is advantageous since the 
voltage—current curve is always monotonous. Therefore, even if 
the reduced model is inaccurate, the gain between the current 
and voltage is always negative. In this case, the target function 
calculates the predicted power by multiplying the voltage and 
current before comparing the result to the desired set point of 
the power: 


J(Au, y) =(1- ett avg Veell = Pset)" W, (avg Veen — Pset) 


$ wer AT WA + Au'SAu (22) 


Again, Jayg and Veen are not calculated in Eq. (22), rather simply 
extracted from Au and y, respectively. 

The nonlinear equality constraint replacing Eq. (20) in this 
case is 


82(Ui, Yi) = Veett,i — FVeey(Ts,i, Mhn,i» -- -) = 0, 
i=0,1...,V=1 (23) 


2.3.3. Inequality constraints 

The changes in the control variables, rather than the actual 
values themselves, are used for the optimization, so a number of 
linear constraints are necessary. First, each control variable has 


a maximum and minimum value. Since the actual variables for 
the optimization problem are defined in terms of changes from 
the previous value, we require, for each time step, that the sum 
of the initial value and the sum of all the preceding steps not 
exceed either the limits for each variable. 

Another constraint places a lower limit on the ratio between 
the inlet fuel flow rate and the current density drawn from the 
fuel cell. Again, since the optimization variables are the changes 
in the control variables this requirement translates into linear 
inequality constraints. For specific details regarding these con- 
straints see Appendix A. 


3. Results and discussion 
3.1. Trade-off between fuel efficiency and performance 


We present results showing the closed-loop response subject 
to a setpoint change from 0.3 to 0.5 W cm7? and back again. 
The control and prediction horizons are five and six time steps, 
respectively, using a time step of 0.5 s. As can be seen, the tem- 
perature of the coolant is damped as is the hydrogen inlet flow 
rate. The output variables are defined as the solid temperature 
and power density. In this example, the temperature is not con- 
trolled at all, so its weight in W, is set to zero. The selection of 
the coefficients of the scaling matrices enable the tuning of the 
overall performance, as will be seen. The MPC tuning param- 
eters used in this study are listed in Table 1. The parameters 
defining the PEM system simulated are identical to those used 
in Golbert and Lewin [11]. 

Fig. 1 shows the closed-loop response obtained when the only 
desire is to ensure robust convergence to the desired setpoint, 
with the value of ®@eff set to zero. As can be seen, the initial 
efficiency is 15%, dropping to 11% as the power approaches 
0.5 Wcm~?, due to the increase in the hydrogen inlet flow rate. 
Subsequently, the current is dropped to lower the power den- 
sity to 0.3 W cm7?, but the fuel flow rate is kept high and the 
efficiency drops to 6.5%. In contrast, Fig. 2 shows the response 
obtained when the efficiency is taken into consideration by set- 
ting werf to 0.2. The first thing to notice is that up to 5 s the system 
is at steady state with an offset in the power density. This is due 
to the trade-off in the optimization target function between con- 
vergence and efficiency (note that the initial efficiency in Fig. 2 
is greater that that of Fig. 1). In the same spirit, the efficiency 


Table 1 

MPC tuning parameters 

Parameter Values 

WwW Since the only output variable specified in 

= the objective function is the tracking error 
from the power set point, the only relevant 
entry scales this error. A value of 5 is used 

w 10!° on the main diagonal 

S The coefficients for the three manipulated 

a variables are 30, 100, and 10, respectively 

At 0.5s 

U 5 

V 6 


J. Golbert, D.R. Lewin / Journal of Power Sources 173 (2007) 298-309 


0 5 10 15 20 2890 
(b) time [s] 


Average 1=0.10339 


5 10 15 20 25 
(c) time [s] 


Fig. 1. Performance with no weight on efficiency and no effective limit on 
cell voltage: (a) Power and voltage trajectories; (b) manipulated variables; (c) 
instantaneous efficiency, 7. 


between 5 and 15s is greater in Fig. 2 at the expense of slower 
convergence to the desired setpoint. This is achieved by using a 
lower fuel flow rate than the previous example. 

Starting from t= 15, however, the fuel flow rate is lowered 
in an attempt to improve the efficiency, although it is obvious 
that the tracking ability of the controller is severely hampered. 


303 


% 5 10 15 20 25290 
(b) time [s] 
Average N=0.11992 
0.16 i r i - 
T 
0.08.5 5 10 15 20 25 
(c) time [s] 


Fig. 2. Performance with weight on efficiency (We = 0.2) and no effective limit 
on cell voltage: (a) Power and voltage trajectories; (b) manipulated variables; 
(c) instantaneous efficiency, 7. 


Close examination of the plot of the current density shows that 
the controller is increasing the current density in an attempt to 
lower the power density—meaning that the controller, either due 
to model inaccuracies or changes in the system conditions, has 
found itself on the right of the peak in the power-—current plot 
(see Fig. 3). 
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P 


g" I 
Fig. 3. Schematic power-current plot. 


In this case, the controller will raise the current density in 
order to lower the power, as required by the change in Pset, 
in Fig. 2. This occurrence has been observed primarily when 
the setpoint is lowered (as opposed to being raised) since both 
considerations (convergence and efficiency) are satisfied by 
lowering the fuel flow rate. Since both criteria tend to lower the 
fuel flow rate, this can lead to excessive changes, which change 
the power—current curve and “trap” the system to the right of 
the peak. When the setpoint is raised, on the other hand, the 
convergence criterion tends to increase the fuel flow rate while 
the efficiency criterion tends to lower it. This conflict of interest 
moderates the change in the fuel flow rate. In an attempt to 
protect the system from crossing the maximum power density, 
we tried defining a lower limit for the voltage and introducing 
this as a boundary for the NLP (remember that the voltage 
itself is one of the optimization variables). The rational for this 
approach can be explained by observing the schematic diagrams 
of the current-power and current-voltage plots in Figs. 3 and 4. 

In essence, the requirement is that the current not exceed the 
value corresponding to the peak of the power plot, /*. Since /* 
corresponds to a specific voltage, limiting the actual voltage to 
be above a limiting voltage, Vmin, will limit the current to the 
left of the power peak. The results using this approach, for the 
case identical to that of Fig. 2, with the cell voltage limited to 
values above 0.5 V, are shown in Fig. 5. 

It is clear that this approach does not solve the problem of 
crossing the power peak. Closer examination of the curves in 
Fig. 6, shown for varying amounts of hydrogen illustrate why 
this is so. Obviously an infinite number of hydrogen/current 
combinations correspond to a given voltage. Thus, setting a limit 
on the voltage implies no constraint on the current, and as seen 
in Fig. 5, the system can indeed jump to the left side of the power 
peak. 

Although it is theoretically possible to define settings that 
could result in better performance, this is not a systematic 
approach and cannot be relied upon as a satisfactory solution. 
A more rigorous solution is to define a nonlinear constraint 


7 
j min 


* 


I I 


Fig. 4. Schematic voltage-current plot. 


w,in 


(b) time [s] 


Fig. 5. Performance with weight on efficiency (We =0.2) and effective limit 
on cell voltage (Veer > 0.5): (a) Power and voltage trajectories; (b) manipulated 
variables. 
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V Increasing 


hydrogen 


I 


Fig. 6. Schematic voltage-current plot for varying hydrogen content. 


whereby the partial derivative of the power to the current is 
above a set value: 


oP 


—>y (24) 
dlayg 


By definition, Eq. (24) ensures that, to the controller’s best 
knowledge, the system is always to the left of the power—current 
peak in Fig. 7. Clearly, if the controller model were perfectly 
accurate, it would be sufficient to ensure that the gradient in 
Eq. (24) is positive (y =0). However, due to model inaccura- 
cies, a larger value for y is used (we have used y =0.2). Fig. 8 
shows the obtained closed-loop response when using the con- 
straint in Eq. (24). Comparing the overall efficiency (12.8%) 
in Fig. 8 to that of Fig. 1 (10.3%) we see that using the effi- 
ciency in the NLP’s target function results in a 25% efficiency 
improvement. 


Increasing 


oo, hydrogen 


P 


I 


Fig. 7. Schematic power-current plot for varying hydrogen content. 


1 
0.9 
0.8 
0.7 


0.6 
0.5 
0.4 


(b) time [s] 


Average n=0.12864 


n n n 1 


5 10 15 20 25 
(c) time [s] 


Fig. 8. Performance with weight on efficiency (We =0.2) and limit on 
OP/dIavyg > 0.2: (a) Power and voltage trajectories; (b) manipulated variables; 
(c) instantaneous efficiency, 7. 


3.2. Ensuring offset-free response 


Although the increase in efficiency is encouraging, the 
offset from the desired power setpoint is troubling. This 
is inherent in the formulation of the NLP, since at steady 
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state, Au=0 and the optimal solution is the minimal value 
of: 


J(Au, y) = (1 — wP — Peet)’ W (P — Pset) 


toa W A (25) 


Is the case of no efficiency consideration weff =0 and the 
optimal solution is complete offset-free convergence. However 
for any other value of wetf, the final value will always display 
offset since the optimum will be some trade-off between the two 
elements of Eq. (25). 

The crux of the problem is the desire to satisfy two — often 
conflicting — criteria, convergence and efficiency. One attempt 
to solve this problem was to add extra weight to the last line 
of the W, matrix. With the intention of forcing the optimizer 
to converge by the end of the prediction horizon. The expecta- 
tion was that the optimizer would find the fuel efficient path to 
convergence by the end of the prediction horizon, if not before 
that. 

However, although the optimizer did propose sequences that 
converged at the end of the prediction horizon, it is important 
to recall that MPC, unlike optimal control, actually imple- 
ments the first step only. Thus, at some point, the optimizer 
decided that, as defined by Eq. (18), it is best to do noth- 
ing at first, reap the benefits of efficiency for a while and 
only later on in the control horizon bring the system to con- 
vergence. If the system is at steady state, and the optimizer 
has decided that no control moves will be made in the fol- 
lowing time step, the subsequent optimizations will be the 
same, ad infinitum, even with offset from the desired set- 
point. Thus, even though the open loop NLP solution predicts 
offset-free convergence, closed loop performance will have off- 
set. 

In searching for a formulation of the NLP that will avoid this 
scenario, we realized that it is crucial to prevent the possibility 
of the optimal solution’s first moves to be zero. This can be 
prevented by formulating a target function that consists only of 
a penalty for control moves, AuT S Au, and imposing a constraint 
on convergence. Since a large number of small control moves 
incur a smaller penalty than a small number of large control 
moves, the optimizer will tend to use all of the moves available 
to it including the first step. With this in mind we defined the 
target function as follows: 


J(Au, y) = wet A(V)? + Au” SAu 26) 


In addition to the equality constraints, Eqs. (19), (21) and, 
(23), we require that: 


P(V) — Pret = 0 en) 


In Eq. (26) only the value of the waste variable, A, at time step V 
(the end of the prediction horizon) appears, and the target func- 
tion is dominated by the penalty on the control moves. Satisfying 
Eq. (27) ensures that the optimizer plans a sequence of control 
moves that converges at the end of the prediction horizon to the 
power setpoint. 
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Fig. 9. Performance with weight on final efficiency, U=5, V=6, (Werf =0.2) 
and limit on 0P/dJzyg > 0.2: (a) Power and voltage trajectories; (b) manipulated 
variables; (c) instantaneous efficiency, n. 


To understand the closed-loop behavior when using this 
formulation, we examine the case where We=0. In this 
case, the solution seeks to minimize the penalty on the 
control moves and results in a monotonous sequence of 
control moves, achieving convergence at the end of the 
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Fig. 10. Performance with weight on final efficiency U=3, V=4, (Werf =0.2) 
and limit on 3P/dIavg > 0.2: (a) Power and voltage trajectories; (b) manipulated 
variables; (c) instantaneous efficiency, n. 


control horizon. It is important to note that for systems 
with more inputs than outputs there are numerous combina- 
tions of control that satisfy the convergence. For example, 
in Fig. 1(b), the system converges to a power setting of 
0.3 W cm~? twice, initially with low hydrogen flow rate and 


subsequently using a higher hydrogen flow rate and a lower 
current. 

For weff =0, any of the possible steady states is satisfac- 
tory. If were is non-zero, however, the system will converge 
to the steady state that satisfies Eq. (27) and has the highest 
efficiency. Unlike before, however, since the target function 
is not affected by the efficiency during the prediction horizon 
the optimal control sequence is monotonous and the dis- 
crepancy between the open and closed loop performance is 
avoided. 

To avoid defining unfeasible constraints, Eq. (27) is combined 
with the objective function of Eq. (26), giving 


J(Au, y) = (PCV) — Peet)” + ep ACV? + AuTSAu (28) 


where ¢ is a very large weight (in the following example, 
y=1000). As p approaches infinity, we approach the original 
formulation of Eqs. (26) and (27). Fig. 9 illustrates the per- 
formance obtained using this new formulation, for the same 
conditions as those in Fig. 8. 

Comparing these results with Fig. 8 clearly shows the 
improved convergence with very little offset (which would 
be completely eliminated if p were infinite). Also note that 
the overall efficiencies for both cases are almost identical. 
Another advantage of this formulation is its ease of tuning. 
Here, the coefficients of S have no affect on convergence to 
the desired setpoint, but only affect the degree to which each 
of the manipulated variables is utilized. Likewise, since the 
convergence is ensured as a constraint, there is no trade-off 
between convergence and efficiency, so tuning We is also 
irrelevant. Instead, the trade-off between the speed of conver- 
gence and efficiency is tuned only by determining the length 
of the prediction horizon, V. This is due to the fact that if 
V is large, there is plenty of time to converge, so the opti- 
mizer can allow itself to aim for an efficient solution. If, 
however, V is shorter, all resources are marshalled to converging 
quickly at the expense of efficiency. Fig. 10 shows the perfor- 
mance of a system identical to that in Fig. 9, with the control 
and prediction horizons reduced to 3 and 4, respectively. As 
expected, the convergence to the desired setpoint exhibited in 
Fig. 10 is faster than that in Fig. 9, at the expense of efficiency 
(11.7% compared with 12.8%). 


4. Conclusions 


As has been previously established, model-based control 
scheme of a PEM fuel cell, relying on a reduced-order, 
nonlinear model of the process, can be used for robust reg- 
ulation. In addition, as demonstrated in this contribution, 
since the controller adjusts a number of manipulated vari- 
ables, it takes advantage of all of the degrees of freedom 
to simultaneously satisfy power demands while optimiz- 
ing the fuel efficiency of the entire system. The use of 
appropriate constraints results in significant improvements 
in fuel efficiency. We have shown that by appropriate for- 
mulation of the objective function, reliable optimization of 
the performance of a PEM fuel cell can be performed in 
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which the main tunable parameter is the prediction and 
control horizons, V and U, respectively. We have demon- 
strated that increased fuel efficiency can be obtained at the 
expense of slower responses, by increasing the values of these 
parameters. 
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Appendix A. Linear inequality constraints 


Each control variable is bounded by a maximum value 
allowed: 


Uj k+1 Umax, 1 
U1 k+U Umax, | 
< : (29) 
Un k+1 Umax,n 
Un,k+U Umax,n 
However 
1 O 0 
U1 k+1 U1 k 
U1 k+U u1,k l 1 0 
= +10 - 0 
Un, k+1 Unk : fey Saas SO 1 0 O 
0 1 
Un, k+U Un,k 
0 0 1 1 
duy x41 Uk 
dui xu Uk 
= + Ay du (30) 
dun k+1 Unk 
dun k+U Unk 
Hence, 


Umax, | ulik 
Umax, 1 ulik 
i | : ee ed |=% (31) 
Umax,n Un,k 
Umax,n Un,k 


In a similar fashion linear constraints are defined for the lower 
boundaries: 


Umin, 1 ] m | 
Umin,1 ulik 
—Aı du < — : —] : =b (32) 
Umin,n Un,k 
Umin,n Un,k 


So far, the constraints ensure that the maximal values of the 
variables will not be exceeded at any step. If the fuel flow rate 
and the current density are to be manipulated there is a danger 
of the optimizer requesting an infeasible current density (above 
the limiting current density, which is largely influenced by con- 
centration overpotential). Thus, for the sake of feasibility (as 
long as the current density is the input to the fuel cell model) a 
minimum ratio between the fuel and the current density must be 
enforced at all times: 


h 
My, 2 9551 = 91 (33) 
where g is a tunable variable, which, in essence, ensures suffi- 


cient saturation of hydrogen. Translating from the values of u to 
the changes in u at each step and defining: 


-1 0 0 0 0 ¢ 0 0 
BS Ne Fe, We. ee ae es a O 
A]! ems he eae, (OY gk Seu ag 
and 
1 
b3 = (UMy,.k — PUL) | : 
1 
gives 
A3 du < b3 (34) 


Note that when defining the matrix, the actual indices depend 
on the number of input variables being used. 
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Similar constraints must be defined for the oxygen/current 
ratio. Combining all of the constraints gives 


Al by 
Al | du < | bo (35) 
A3 b3 


This defines all of the constraints. For the sake of sensitivity, 
the variables are all scaled by their respective values entering 
the optimization: 


du; = du; ui nom (36) 
or 
U1,nom 
U1,nom 
du=diag} : | du* (37) 
Un, nom 
Un, nom 


Substituting Eq. (37) into Eq. (35) gives 


U1 nom 
AL U1 nom by 
=A - diag du* < | b2 (38) 
A3 Un,nom b3 


Un,nom 


These are the linear constraints on the optimization vari- 
ables. Furthermore, in cases where the hydrogen or oxygen flow 
rates are not optimized, the current is limited to the permitted 
ratio between the current and the constant value of the reactant 
flow rate. In this case, the upper limit on the current is either 
the set maximum current density (imposed by the user) or the 
value determined by the permitted current/reactant flow ratio, 
the lower of the two. 

Note that since the optimization variables include the state 
and output variables as well, the constraint matrix needs to be 
padded with zeros to match the dimensions of the optimization 
variable. 
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